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SUMMARY  PAGE 


THE  PROBLEM 

We  are  using  event-related  potentials  (ERPs)  to  characterize  the  changes 
in  cognition  that  occur  when  a  high  level  of  performance  must  be  sustained  for 
many  hours  despite  fatigue  and  loss  of  sleep.  The  need  for  such  performance 
occurs  during  the  repeated,  long-range  tactical  missions  that  characterize 
sustained  and  continous  operations  in  naval  aviation. 

FINDINGS 

We  describe  a  computer  program  we  have  developed  that  scans  ERP  data  for 
several  types  of  electrical  artifacts.  When  feasible,  the  program  attempts  to 
correct  the  electroencephalogram  (EEG)  for  the  effects  of  eye  movements  and 
blinks.  After  doing  so,  the  program  checks  the  success  of  its  corrections, 
writes  a  data  file  of  the  remaining  artifacts,  and  writes  the  corrected  EEG  to  a 
data  file.  A  copy  of  the  program's  source  code  is  included  as  an  appendix. 

RECOMMENDATIONS 

Clearly,  ERP  data  must  be  routinely  screened  for  electrical  artifacts. 

The  algorithms  used  in  this  program  represent  several  contemporary  approaches 
to  the  artifact-detection  problem.  None  of  these  algorithms  is  fully 
satisfactory  when  used  in  isolation.  Combining  them,  however,  yields  an 
improved  quality-control  system  for  ERP  studies. 
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INTRODUCTION 


An  ERP  is  an  aggregate  electrical  field  produced  when  a  population  of 
neurons  discharge  simultaneously  in  response  to  a  sensory  or  internal  event. 
The  amplitudes  of  ERPs  are  greatly  attenuated  as  they  propagate  from  the  brain 
to  the  surface  of  the  scalp;.  As  a  result,  ERPs  recorded  at  the  scalp  are 
easily  contaminated  by  stray  electrical  events.  These  electrical  artifacts 
can  render  data  sets  unusable.  If  they  pass  undetected,  they  can  render  a 
data  set  misleading.  Consequently,  artifact  detection  is  among  the  most 
important  aspects  of  ERP  recording  technique.  For  these  reasons,  we  have 
written  a  computer  program  that  scans  ERP  data  for  several  types  of  electrical 
artifacts.  This  report  describes  our  program. 

We  review  the  topic  of  electroencephalogram  ''EEG)  and  ERP  artifacts  only 
briefly;  a  comprehensive  discussion  of  the  topic  can  be  found  elsewhere  (1). 

We  use  the  term  artifact  here  to  refer  to  any  potential  that  makes  a  sample  of 
EP.P  data  unusable,  including  random  noise. 

The  amplitudes  of  ERPs  recorded  with  scalp-surface  electrodes  are  small, 
within  roughly  an  order  of  magnitude  of  1  uV.  Hence,  they  are  easily 
contaminated  by  extraneous  voltages.  Such  artifacts  may  be  biological 
potentials  of  nonneuronal  origin  or  nonbiological  potentials  generated  by  the 
recording  system  or  nearby  electrical  equipment.  For  example,  artifacts  in 
ERP  recordings  include  potentials  from  eye  movements  and  blinks,  muscle 
contractions,  tongue  movements  (the  tongue  is  polarized  end-to-end),  galvanic 
skin  responses,  changes  in  the  skin-electrode  interface,  and  movements  of 
electrode  leads.  Line-power  artifacts  occur  when  equipment  is  improperly 
grounded,  electrode  impedances  are  high,  or  electrode  leads  are  placed  too 
near  power  cords  or  other  radiant  sources.  Aliasing  artifacts  occur  when  the 
EEG  is  converted  to  a  digital  format  without  first  removing  frequency 
components  above  one-half  of  the  sampling  rate  (1). 

Adjusting  instrumentation  and/or  procedures  will  minimize  or  eliminate 
many  artifacts.  Others  are  difficult  or  impossible  to  avoid.  Artifacts  from 
subjects  are  usually  least  controllable:  Subjects  blink,  move,  and  emit 
galvanic  skin  responses  despite  instructions  to  the  contrary.  Usually,  the 
best  one  can  do  is  to  identify  any  artifacts  and  ensure  that  they  do  not 
compromise  the  data.  Our  computer  program  should  assist  in  that  effort. 

This  report  contains  three  major  sections  ard  an  appendix.  Section  1 
discusses  artifact  detection  and  correction  techniques  of  the  artifact 
filtering  program,  Artfil;  section  2  describes  how  the  program  is  used;  ana 
section  3  contains  system  requirements.  The  appendix  includes  a  copy  of  the 
program's  source  code. 

1.  ARTIFACT  DETECTION  AND  CORRECTION 


ARTIFACT  DETECTION 

Artfil  checks  the  EEG  and  electrooculcgram  (EOG)  for  six  types  of 
artifacts:  (1)  eyeblinks,  (2)  voltage  spikes,  (3)  large  absolute  voltages  (4) 

large  root-mean-square  (rms)  voltages,  (5)  dead  EEG  channels,  and  (6) 
amplifier  clipping  (or  saturation) . 
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1.  Ere  blink  detection  is  performed  using  an  algorithm  described  by 
Gratton,  Coles,  and  Donchin  (2).  A  blink  is  identified  when  the  local  slope 
of  the  vertical  EOG  trace  exceeds  a  critical  value,  blinkcrit:.  A  local  EOG 
slope  exceeding  blinkcrit  indicates  eyelid  movement  and,  hence,  a  blink. 

After  Gratton  et  al . ,  we  define  blinkcrit  as  a  criterion  change  in  voltage 
during  a  10-ms  interval. 

2.  Spike  detection  is  performed  similarly  to  blink  detection.  Spikes 
are  identified  when  the  difference  between  successive  voltage-values  in  an  EEC 
or  EOG  time  series  exceeds  a  criterion  value,  spikecrit .  The  vertical  EOG  is 
not  scanned  for  spikes,  which  could  be  confused  with  blinks.  If  an  EEG 
channel  or  the  lateral  EOG  channel  contains  a  spike  artifact,  the  data  are  not 
corrected  for  ocular  artifacts  and  are  flagged  for  rejection. 

3.  Large  absolute  voltages  are  detected  by  comparing  each  voltage  in 
the  EEG  time  ser’es  to  a  criterion  voltage,  artcr it .  after  the  eye-movement 
compensation  algorithm  described  presently  has  been  applied  to  the  EEG  data. 
Any  epoch  and  channel  of  EEG  data  that  contains  an  absolute  voltage  exceeding 
artcri t .  after  ocular  artifact  compensation,  will  be  flagged  for  rejection. 
This  is  probably  the  most  widely  used  technique  for  detecting  artifacts  in  ERP 
data.  A  check  of  this  type  is  often  applied  to  a  single  EOG  channel  or  a 
frontal  EEG  channel.  In  contrast,  Artfil  checks  all  EEG  channels  for  absolute 
voltage  artifacts  after  correcting  the  EEG  for  EOG  contamination.  It  does  not 
search  the  EOG  for  absolute -voltage  artifacts  of  this  type.  It  does,  however, 
search  the  EOG  for  voltages  large  enough  to  cause  amplifier  saturation  because 
amplifier  saturation  renders  ocular-artifact  compensation  impossible. 

4.  Large  rms  voltages  are  detected  by  comparing  the  overall  root- 
mean-sq aared  amplitude  of  each  epoch  of  data  with  a  criterion  value,  rmscrit . 
The  quantity  compared  to  rmscrit  is  the  standard  deviation  of  the  points 
within  a  given  epoch  and  channel.  An  epoch  and  channel  of  data  will  be 
flagged  for  rejection  if  its  rms  voltage  exceeds  rmscrit  after  ocular  artifact 
compensation  has  been  performed.  The  rmscrit  should  be  set  to  a  value  sub¬ 
stantially  smaller  than  the  absolute -voltage  criterion,  voltcri t .  This  is 
because  overall  rms  voltages  are  inherently  smaller  than  peak  voltages. 

5.  Dead  EEG  channels  are  detected  by  comparing  the  overall  rms 
amplitude  of  each  epoch  of  data  with  a  criterion  value,  deader it .  Any  epoch 
and  channel  with  an  overall  rms  voltage  less  than  deadcrit  will  be  flagged  for 
rej  ection. 

6.  Amplifier  saturation  is  detected  by  comparing  the  absolute  value 
of  each  voltage  in  each  epoch  of  data  with  a  criteric  -  lue ,  cllpcrlt .  For 
computational  speed  and  convenience,  we  define  cliucrit  as  an  input  voltage 
large  enough  to  cause  an  amplifier  to  saturate.  That  is,  the  output  voltage 
that  will  produce  clipping  is  divided  by  amplifier  gain  to  yield  the  value  of 
clipcrit.  A  error  in  amplifier  calibration  could  cause  Artfil  to  miss 
instances  of  clipping.  This  can  be  minimized  by  setting  clipcrit  to  90-95%  of 
an  amplifier's  actual  clipping  voltage,  assuming  calibration  accuracy. 

Separate  clipping  criteria  can  be  applied  to  the  EEG  and  EOG  data,  as 
discussed  in  section  2. 
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TREATMENT  OF  OCULAR  ARTIFACTS 


Eye  movements  and  blinks  produce  electrical  potentials  that  can  be 
recorded  with  EEG  electrodes.  Changes  in  the  spatial  distribution  of  the 
eyes'  standing  electrical  fields  that  occur  when  the  eyes  move  in  the  head 
cause  eye-movement  artifacts.  Eyeblink  artifacts  are  transient  changes  in 
these  fields  due  to  resistance  changes  associated  with  eyelid  movements  (see 
reference  1  for  a  discussion) .  The  locations  of  the  active  and  reference 
electrodes  determine  the  amplitudes  of  both  types  of  artifacts.  When  recorded 
by  EEG  electrodes  placed  at  sites  commonly  used  for  ERP  recording,  the  ampli¬ 
tudes  of  ocular  artifacts  are  often  larger  than  those  of  ERPs. 

Generally,  eye  movements  are  monitored  by  recording  at  least  one  channel 
of  EOG  along  with  the  EEG  to  detect  EEG  segments  that  may  be  contaminated  by 
ocular  potentials.  The  customary  procedure  is  to  reject  an  epoch  of  EEG  when 
the  EOG  exceeds  a  criterion  absolute  voltage.  This  widely  used  strategy  is 
simple  and  requires  only  two  assumptions.  The  first  assumption  is  that  the 
detection  procedure  is  sensitive  enough  that  any  contamination  produced  by 
undetected  eye  movements  can  be  safely  ignored.  The  second  is  that  excluding 
EEG  with  eye-movement  artifacts  does  not  produce  a  biased  data  set.  To  our 
knowledge,  neither  assumption  has  been  thoroughly  examined  and  verified. 

An  alternative  is  to  use  EOG  recordings  to  remove  ocular  potentials  from 
the  EEG  (see  reference  1  for  a  review),  as  in  Artfil.  This  approach  also 
involves  assumptions,  which  center  on  the  accuracy  with  which  direct 
recordings  of  the  EOG  can  be  used  to  estimate  the  EOG  contamination  present  in 
EEG  recordings.  All  compensation  procedures  of  this  type  involve  subtracting 
suitably  scaled  EOG  waveforms  from  the  EEG  waveforms.  (Sometimes  the  waves 
are  decomposed  into  frequency  components,  and  the  different  frequencies  are 
scaled  separately.)  Because  the  EOG  is  recorded  from  the  head,  EOG  electrodes 
probably  always  record  some  EEG  activity.  Hence,  these  procedures  probably 
distort  ERP  data  somewhat  because  they  involve  subtracting  brain  activity 
recorded  by  the  EOG  electrodes  from  brain  activity  recorded  by  the  EEG 
electrodes.  The  magnitude  of  this  problem  has  not  been  thoroughly  studied. 

One  way  to  reduce  EOG  contamination  is  to  record  the  EOG  differentially  from  a 
pair  (or  pairs)  of  electrodes  adjacent  enough  so  ;hat  the  local  EEG  is  nearly 
identical  in  each  (3).  Differential  amplification  will  then  tend  to  remove 
the  EEG  from  the  EOG  recording. 

The  ocular  artifact  compensating  routines  of  Artfil  presume  that  bipolar 
EOG  recordings  result  from  two  pairs  of  electrodes,  each  containing  a  record¬ 
ing  electrode  and  a  reference  electrode.  The  program  assumes  that  one  pair  of 
electrodes  obtained  vertical  EOG  data,  perhaps  from  an  electrode  above  one  eye 
referred  to  an  electrode  below  that  eye.  These  data  are  used  to  detect  and 
compensate  the  data  for  eyeblink  artifacts.  The  program  also  assumes  a  second 
pair  of  electrodes  obtained  horizontal  or  oblique  EOG  data.  We  refer  to  this 
electrode  pair  as  the  lateral  EOG  channel.  Such  data  might  be  recorded  from 
electrodes  placed  to  the  left  and  right  of  one  eye,  or  obliquely  about  one 
eye.  Artfil  uses  these  data  to  correct  the  EEG  for  eye  movements. 

Artifil  corrects  the  EEG  for  artifacts  caused  by  eye  movements  and  blinks 
using  a  variant  of  the  procedure  described  by  Gratton  et  al .  (2).  Their 
algorithm  assumes  that  the  waveform  recorded  from  an  EEG  electrode  can  be 
approximated  as  the  sum  of  two  time  series.  The  first  time  series  is  the 
actual  EEG  waveform;  the  other  is  a  linearly  attenuated  version  of  the  EOG. 


By  this  assumption,  the  actual  EEG  can  be  recovered  from  the  recorded  EEG  by 
subtracting  the  appropriately  scaled  EOG  point-by-point . 

The  basic  procedure  consists  of  two  steps,  each  carried  out  on  an  epoch- 
by-epoch  basis,  In  the  first  step,  a  least-squares  cross  regression  is 
calculated  using  the  EEG  time  series  as  criterion  variables  and  the  EOG  as 
predictor  variables.  This  regression  estimates  the  constant  to  multiply  the 
EOG  time  series  to  get  the  best  linear  prediction  of  the  EEG  time  series.  The 
constant  is  an  estimate  of  the  proportion  of  the  EOG  contained  in  the  recorded 
EEG.  Its  value  is  estimated  separately  for  each  EEG  recording  site.  In  the 
second  step,  part  of  the  recorded  EOG  is  subtracted  in  pointwise  fashion  from 
the  EEG.  The  amount  subtracted  from  the  EEG  is  determined  by  the  proportiona¬ 
lity  constant  estimated  in  step  one. 

These  procedures  are  different  for  eye  movements  and  blinks.  For  eye 
movements,  the  proportion  is  estimated  from  complete  EEG  and  EOG  time  series. 
For  blinks,  the  proportion  of  EOG  subtracted  from  the  EEG  is  estimated  from 
data  obtained  during  periods  in  which  the  local  slope  of  the  vertical  EOG 
(calculated  in  a  10-ms  time  window)  exceeds  a  criterion  value. 

When  the  blink-slope  criterion  is  properly  chosen,  the  data  used  to 
estimate  the  proportionality  constant  for  eyeblink  compensation  will  be 
selected  from  the  leading  and  trailing  edges  of  eyeblink-artifact  waveforms 
(i.e.,  when  the  eyelids  are  moving  rapidly).  Selecting  an  appropriate  value 
for  this  criterion  is  critical  to  the  performance  of  the  algorithm.  If  the 
blink-scope  criterion  is  too  high,  the  algorithm  either  fails  to  detect  blinks 
or  produces  unstable  estimates  of  the  blink-scaling  constant  based  on  small 
numbers  of  data  points.  On  the  other  hand,  if  the  value  is  too  low,  the 
algorithm  mistakes  EOG  noise  for  blinks  and,  thereby,  underestimates  the 
blink-scaling  constant  and  undercorrects  the  EEG  for  the  effects  of  blinks. 

We  can  only  advise  on  how  to  select  a  value  for  the  blink  criterion.  The 
shapes  of  eyeblink  waveforms  depend  on  where  the  EOG  electrodes  are  placed  and 
on  how  the  data  are  filtered.  Blinks  also  vary  from  one  individual  to  another 
and  from  one  blink  to  the  next.  We  set  the  slope  criterion  to  a  value  that  is 
less  than  the  slopes  of  the  blink  waveform  leading  edg>s  when  blinks  are  25- 
75%  of  their  maximum  amplitudes.  Because  blinks  are  variable  and  because  the 
results  of  different  criterion  settings  must  be  checked  visually  blink  by 
blink,  selecting  a  value  of  the  blink  criterion  can  be  time-consuming. 

The  procedure  to  compensate  for  ocular  artifacts  described  by  Gretton  et 
al.  (2)  includes  a  third  step  that  is  not  included  in  Artfil.  It  involves 
subtracting  a  signal-average  of  all  epochs  of  the  ERP  data  (with  ocular 
artifacts  included)  from  each  individual  epoch  of  data  before  estimating  the 
scaling  constants  and  subtracting  the  scaled  EOG  from  the  EEG.  This  step  is 
intended  to  remove  the  brain  activity  evoked  by  eyeblinks  from  the  EEG  data. 

We  examined  the  procedure  and  found  that  in  our  hands  the  ocular-artifact 
compensating  algorithm  performed  better  without  it. 
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2.  USING  THE  PROGRAM 


DATA  FILES 

Artfil  requires  two  input  data  files  and  creates  two  output  data  files. 

The  main  input  data  file  contains  digitized  EEG  and  EOG;  the  default  name  is 
■  pochchan.  The  second  input  data  file  contains  parameters  used  by  Artfil  for 
artifact  detection  and  EOG  artifact  compensation.  The  default  name  of  this 
file  is  artifact. p2. 

The  main  output  data  file,  the  default  name  of  which  is  epochchan.c, 
contains  digitized  EEG  with  the  ocular-artifact  compensating  algorithm 
applied.  Because  Artfil  does  not  discard  data  with  uncorrectable  artifacts, 
epochchan.c  contains  as  much  EEG  data  as  the  input  file  epochchan.  Artifil 
writes  the  epoch  and  channel  numbers  of  data  with  uncorrectable  artifacts  to  a 
second  output  data  file  named  artifacts. 

The  physiological  data  in  epochchan  and  epochchan.c  are  stored  as  time 
series  of  digitized  EEG  and  EOG  amplitudes  scaled  in  0.1-uV  units.  Values  of 
the  EEG  and  EOG  data  are  stored  in  the  files,  which  are  direct-acce  .s,  as  two- 
byte  integers.  Each  time  series  corresponding  to  data  from  one  recording 
channel  in  one  epoch  is  stored  as  a  separate  record. 

The  layout  of  the  epochchan  files  is  best  illustrated  by  considering  each 
voltage  in  the  file  as  an  element  in  a  triply  subscripted  voltage  array, 
Y(e.£.,i).  The  subscript  e  indexes  the  ordinal  (and  temporal)  position  of  the 
recording  epoch  from  which  y  was  obtained.  The  value  of  e  varies  from  1  to 
ne .  the  number  of  epochs  in  the  data  set  under  consideration.  The  index  c 
indicates  the  number  of  the  recording  channel  from  which  v  was  obtained.  The 
value  of  c  varies  from  1  to  n£,  the  number  of  recording  channels.  The  index  i 
refers  to  the  ordinal  (and  temporal)  position  of  v  within  the  current  record¬ 
ing  epoch.  The  value  of  i  varies  from  1  to  ni.  Thus,  v(l,2,3)  is  the  third 
voltage  point  obtained  from  channel  two  in  epoch  one. 

Values  of  v  are  organized  in  the  epochchan  files  such  that  index  j.  varies 
most  rapidly.  (Recall  that  Qi  voltages  from  channel  c  in  epoch  e  comprise  one 
record  of  epochchan.)  The  channel  index,  c,  varies  next  most  rapidly,  and  the 
epoch  index,  e,  varies  least  rapidly. 

The  input  data  file  artifact. p2  contains  several  parameters  that  Artfil 
uses  to  detect  artifacts.  The  user  should  tailor  the  parameters  for  each 
recording  system  and  experiment  and  enter  one  number  per  line  in  the  file  as 
ASCII-coded  numbers  in  the  order  indicated  below. 

clipcrit  -  The  criterion  voltage  for  detecting  EEG  amplifier 
clipping.  For  simplicity,  the  value  of  clipcrit  is  expressed  as  the  voltage, 
in  uV  at  the  amplifier  inputs,  sufficient  to  produce  clipping.  Any  epoch  of 
data  in  which  clipping  is  detected  will  not  be  corrected  for  ocular  artifacts 
ant  will  be  flagged  for  rejection.  A  conservative  value  of  clipcrit,  assuming 
EOG  amplifier  gains  of  20,000  and  saturation  voltages  of  5.0  V,  would  be 
225.0.  A  separate  voltage  criterion  is  used  to  detect  EOG  amplifier  clipping 
(see  clipcrit2  below). 
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cliocrit2  -  A  second  clipping- criterion  voltage  that  is  used  to 
detect  EOG  amplifier  clipping.  Again,  Artfil  detects  amplifier  clipping  by 
monitoring  amplitudes  expressed  as  amplifier- input  voltages.  The  amplifier- 
input  voltages  of  the  EOG,  however,  are  much  larger  than  those  of  the  EEG. 
Hence,  the  EOG  and  EEG  amplifiers  may  be  operated  at  different  gains.  If  the 
EOG  and  EEG  amplifiers  differ  in  gain  settings  (but  are  otherwise  similar) 
they  will  saturate  at  different  input  voltages. 

A  solution  to  this  problem  is  to  define  separate  clipping  criteria 
for  the  EEG  and  EOG  amplifiers:  Artfil  uses  the  value  of  clipcrit2  to  detect 
EOG  amplifier  clipping.  Any  epoch  in  which  the  absolute  voltage  of  the  EOG 
exceeds  clipcrit2  will  not  be  used  in  ocular-artifact  correction  and  will  be 
flagged  for  rejection.  A  conservative  value  of  clipcrit2,  assuming  EOG  ampli¬ 
fier  gains  of  2000  and  output  saturation  voltages  of  5.0  V,  would  be  2250.0. 

blinkcrit  -  A  criterion  for  detecting  blinks.  After  the  method  of 
Gratton  et  al.  (2),  blinkcrit  is  defined  as  a  criterion  value  of  the  local 
slope  of  the  vertical  EOG  trace  measured  in  a  10-ms  time  interval.  Artfil 
adjusts  the  input  value  of  blinkcrit  linearly  to  accommodate  sampling  inter¬ 
vals  that  differ  from  10  ms.  As  discussed  previously,  determining  an  appro¬ 
priate  value  for  blinkcrit  may  require  some  experimentation. 

solkecrit  -  A  rate  of  voltage  change,  in  uV/ms ,  that  identifies  the 
presence  of  a  high-frequency  spike  in  an  EEG  channel.  The  criterion  is 
applied  to  differences  calculated  between  successive  voltages  in  each  epoch 
and  channel  of  data  except  for  the  vertical  EOG.  As  noted,  the  EOG  is  not 
scanned  for  spikes  for  fear  of  confusing  them  viith  blinks.  A  suggested  value 
of  spikecrit  is  33.0. 

voltcrit  -  An  absolute-voltage  criterion,  expressed  in  uV,  used  for 
rejecting  epochs  of  data.  Any  epo-.h  and  channel  of  EEG  data  that  contains  an 
absolute  voltage  exceeding  voltcrit  after  the  ocular  artifact  compensation 
procedure  has  been  applied  will  be  flagged  for  rejection.  A  suggested  value 
of  voltcrit  is  40  uV. 

rmscrit  -  A  root-mean-squared  voltage  criterion  used  for  rejecting 
epochs  of  data  (scaled  in  uV) .  The  standard  deviation  of  the  points  within  a 
given  epoch  and  channel  are  compared  to  rmscrit .  An  epoch  and  channel  of 
data  will  be  flagged  for  rejection  if  its  rms  voltage  exceeds  rmscrit  after 
ocular  erti Tact  compensation  has  been  performed. 

deader i u  -  A  root-mean- squared  voltage  criterion  used  to  identify 
dead  amplifiers  (scaled  in  uV) .  Any  epoch  and  channel  with  an  overall  rms 
voltage  lejs  than  dsadcrit  will  be  flagged  for  rejection.  We  suggest  a  value 
near  4.0  uV. 

samprate  -  The  EEG  and  EOG  sampling  rate  scaled  in  data  points  per 
second.  Artfil  uses  samprate  when  applying  the  EOG  slope  criterion 
(blinkcrit)  to  the  data.  This  variable  is  an  integer. 

The  output  file  artifacts  is  a  summary  of  the  artifacts  found  by  Artfil 's 
artifact  detecting  algorithms.  The  file  is  written  to  the  directory  in  which 
the  program  found  the  original  epochchan  file.  the  file  contains  one  record 
of  artifact  data  for  each  epoch  and  channel  of  data  examined.  As  in  epochchan 
files,  data  records  for  different  channels  are  written  in  ordered  blocks. 
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Each  block  of  records  contains  the  data  from  channels  for  a  single  epoch. 
Each  record  within  a  block  corresponds  to  the  data  from  one  channel. 

Each  analytical  program  uses  the  information  in  artifacts  differently. 
Generally,  a  program  opens  epochchan.c  and  artifacts  files  simultaneously.  It 
then  reads  the  information  in  the  kth  record  of  Artifacts,  which  is  the 
summary  of  artifacts  for  the  physiological  data  stored  in  record  k  of 
epochchan.c.  The  program  uses  this  information  to  determine  whether  the  epoch 
of  data  in  record  k  of  epochchan.c  is  usable.  If  so,  the  program  reads  the 
physiological  data  and  analyzes  it.  If  not,  the  program  either  proceeds  to 
another  record  or  exits  as  appropriate.  Each  record  in  artifacts  contains 
ten,  2-byte-integer  data  fields  described  below. 


Field  Interpretation 

1.  1  if  a  clip  or  spike  was  detected  in  the  EOG  in  the  current  epoch,  or 
if  a  clip,  spike,  rms ,  or  absolute -voltage  artifact  was  detected  in 
the  EEG  in  the  current  epoch;  0  otherwise. 

2.  The  number  of  the  current  epoch. 

3.  The  number  of  the  current  channel. 


4. 


5. 


6. 


7. 


8. 


1  if  the  current  channel  is  the  vertical  EOG  and  a  blink  was  detected 
in  the  current  epoch;  0  otherwise. 


1  if  the  current  channel  is  an  EEG  channel  and  an  absolute -voltage 
artifact  was  detected  in  the  current  epoch;  0  otherwise. 


1  if  the  current  channel  is  an  EEG  channel  and  an  rms  artifact  was 
detected  in  the  current  epoch;  0  otherwise. 

1  if  the  current  channel  is  dead  in  the  current  epoch;  0  otherwise. 


1  if  the  current  channel  is  not  the  vertical  EOG  channel  and  a  spike 
artifact  was  found  in  the  current  epoch;  0  otherwise. 


9.  1  if  a  clipping  artifact  was  found  in  the  current  channel  and  epoch; 

0  otherwise. 


10.  The  overall  rms  amplitude  of  the  current  channel  and  epoch,  rounded 

to  the  nearest  uV. 


RUNNING  THE  PROGRAM 

Artfil  allows  the  user  to  specify  several  arguments  on  the  command  line 
when  the  program  is  loaded.  To  run  Artfil,  the  user  enters: 

programname  ARGUMENTS 

where  programname  is  the  name  assigned  to  the  executable  version  of  the 
program,  and  ARGUMENTS  is  a  sequence  cf  command  line  arguments.  Only  the 
first  two  arguments  (-c  and  -n)  are  required;  the  rest  are  optional.  The  list 
of  valid  command  line  arguments  is: 
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-cCnumber  of  channels>,  required 

The  number  of  channels  of  data  in  the  epochchan  file. 

-n<number  of  points  per  epoch>,  required 

The  number  of  data  points  per  channel  per  epoch. 

[-vCvertical  ey-*  channel>]  ,  default  -  1 

The  channel  number  of  the  vertical  EOG  channel  if  not  channel  1. 

[-Klateral  eye  channel>]  ,  default  -  2 

The  channel  number  of  the  lateral  EOG  channel  if  not  channel  2. 

[-F<first  epoch  to  process>] ,  default  -  1 

The  first  epoch  processed  in  the  current  run.  If  too  large  to 
process  at  one  time,  an  epochchan  file  can  be  processed  as 
blocks  of  epochs  in  sequential  runs.  The  first  and  last  epoch 
of  each  run  indicated  with  this  argument  and  the  -L  argument. 

If  not  specified,  the  first  epoch  processed  will  be  epoch  1. 

[-L<last  epoch  to  process>] ,  default  -  10 

The  last  epoch  processed  in  the  current  run.  This  argument  is 
used  in  conjunction  with  the  -F  argument.  If  not  specified,  the 
last  epoch  processed  will  be  epoch  10. 

[-e(chan.  to  omit  from  artifact  filtering)],  default  -  no  omissions 
Indicates  which  channel  of  data  should  not  be  processed  by  the 
artifact  filter.  Data  in  any  channel  so  indicated  will  be 
written  to  the  output  file  epochchan. c  exactly  as  it  was  read 
from  the  input  file  epochchan.  The  argument  may  be  listed  more 
than  once  on  a  given  command  line. 

[-m  (if  present,  perform  three-point  filtering  on  data)],  default  - 

none 

Used  for  three -point  smoothing  of  the  input  data. 

[ -V(verbose) ] ,  default  -  no  verbose 

If  this  argument  is  present,  Artfil  will  print  numerous  messages 
pertinent  to  the  status  of  the  analysis. 

[-P<artifact  parameter  file  path  name>] ,  default  -  artifact, p2 

Used  if  the  artifact  parameter  file  is  not  the  default.  Enter 
the  new  file  name. 

[-Icinput  file  path  name>] ,  default  -  epochchan 

Used  when  the  input  physiological  data  file  is  not  the  default. 
Enter  the  new  file  name. 

(-OCoutput  file  path  name>] ,  default  -  epochchan, c 

Used  when  the  output  physiological  data  file  is  not  the  default. 
Enter  the  new  file  name. 

[-A<artifact  output  file  path  name>],  default  -  artifacts 

Used  when  the  artifact  summary  file  is  not  t1  default. 

Enter  the  new  file  name. 
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For  example,  entering  'artfil'  will  display  a  help  screen  containing  a 
list  of  command  line  options.  If  we  then  enter: 

artfil  -clO  -n250  -vl  -12  -F125  -L250  -elO 

the  program  will  read  10  channels  of  data  (-clO)  comprising  250  data  points 
'-er  channel  per  epoch  (-n250).  The  program  will  treat  channel  1  as  the 
vertical  EOG  (-vl)  and  channel  2  as  the  lateral  EOG  (-12).  Artfil  will 
analyse  epoch  125  (-F125)  first  and  epoch  250  (-L250)  last  but  will  not 
examine  data  in  channel  10  (-r.lo/.  If  no  data  files  are  specified  on  the 
command  line,  the  default  files  will  be  used  by  the  program. 


3.  SYSTEM  REQUIREMENTS 


DATA  LIMITATIONS 

For  efficiency  on  the  MASSCOMP  5500,  the  maximum  number  of  channels  is 
20,  the  maximum  number  of  epochs  per  run  is  350,  and  the  maximum  number  of 
points  per  channel  per  epoch  is  400.  These  values  were  set  to  allow  Artfil  to 
run  entirely  with  nonvirtual  arrays.  They  can  be  modified  depending  on 
available  memory.  To  adjust  the  maximum  channels,  epochs,  or  points,  the 
following  parameter  statements  can  be  edited  in  Artfil 's  source  code: 

parameter  (MAXCHANNELS  -  20)  !  Maximum  number  of  channels 

parameter  (MAXEPOCHS  -  350)  !  Maximum  number  of  epochs 

parameter  (MAXPOINTS  -  400;  !  Maximum  number  of  points  per 

epoch 

COMPATIBILITY 

Artfil  runs  on  a  MASSCOMP  5500  computer  under  MASSCOMP  Real-Time  Unix 
Version  4.0.  It  is  written  in  Fortran  77  with  DEC  VAX  extensions,  and  it 
should  be  reasonably  compatible  with  most  compilers  on  VAX,  Hewlett-Packard, 
IBM,  and  IBM  PC-ccmpatible  computers, 

The  major  departure  from  standard  Fortran  77  in  Artfil  is  the  use  of  do- 
loops  that  end  with  "end do"  statements.  Some  compilers  do  not  support  control 
structures  of  this  type.  They  can  be  replaced  easily  with  traditional  loops. 
Another  difference  from  standard  Fortran  77  is  the  use  of  command  line  argu¬ 
ments  in  Artfil.  If  a  compiler  does  not  support  command  line  arguments, 

Artfil  can  be  modified  to  read  the  argument  values  from  a  file. 

Lastly,  before  the  variable  declaration  statements,  Artifil  uses  an 
"implicit  none"  statement,  which  requires  that  all  variables  be  explicitly 
declared.  It  can  be  removed.  Other  incompatibilities  may  occur  in  input  and 
output  statements,  which  are  notoriously  nonstandard  across  Fortran  compilers. 
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APPENDIX 


program  artfil 
c 

c  Author:  Dr.  Robert  R.  Stanny  (NAMRL) ,  September  1986 

c 

c  Modification  Log: 
c 

c  1.  November  1987,  Sam  J.  La  Cour,  Jr.  (NARDAC) 
c 

c  Made  the  program  compatible  with  the  files  and  methods  used  in 

c  data  acquisition  programs  written  for  the  system  circa  1987. 

c 

c  2.  January  1989,  Sam  J.  La  Cour,  Jr.  (NARDAC) 

c 

c  Involved  allowing  specification  of  eye  channels 

c  and  changing  the  way  in  which  blinks  are  scaled 

c  prior  to  calculation  of  scaling  constant, 

c 

c  3.  September  1989,  Sam  J.  La  Cour,  Jr.  (NARDAC) 
c 

c  Place  "clipcrit2"  variable  in  parameter  file  after  "clipcrit". 

c 

c  4.  January  1990,  Sam  J.  La  Cour,  Jr.  (NARDAC) 

c 

c  Modified  calculation  of  adjusted  blink  and  spike  criteria, 

c 

c  End  Modifications 
c 

c  Program  Usage : 
c 

c  1 .  Command  Line : 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


artfil  -cCnumber  of  channels>,  required 

( -vCvertical  eye  channel>] ,  default  1 
(-Klateral  eye  channel>)  ,  default  2 
[-F<first  epoch  to  process>] ,  default  1 
[-LClast  epoch  to  process>] ,  default  10 
[-eCchannel  to  omit  from  artifact  filtering>] , 
default  no  omissions 

[ -m  (if  present,  perform  3  point  filtering  on  data)], 
default  none 

[ -V  (verbose)],  default  no  verbose 
[-P<artifact  parameter  file  path  name>] , 
default  "artifact , p2" 

[-I<input  file  path  nar,e>]  ,  default  "epochchan" 
[-OCoutput  file  path  name>] ,  default  "epochchan . c" 

[ -A<artifact  output  file  path  name>] , 
default  "artifacts" 


c 

c  2.  The  program  assumes  the  existance  of  "epochchan"  in  the  current 
c  directory.  This  file  contains  the  EOG  and  EEG  data.  Also,  the 
c  file  "epochinfo"  must  exist  in  the  same  directory.  It  contains 
c  the  number  of  eppchs  (first  line)  and  the  number  of  points  per 
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c  epoch  (second  line).  This  is  a  text  file, 
c 

c  3.  The  program  reads  the  following  parameters  from  the  artifact 
c  parameter  file: 
c 

c  clipcrit  uV  voltage  used  to  detect  amplifier  clipping  in  EEG  channels 

c  clipcrit2  uV  voltage  used  to  detect  amplifier  clipping  in  EOG  channels 

c  blinkcrit  uV/lOms  voltage  change  used  as  a  criterion  for 

c  identifying  blinks  in  eye  channels 

c  spikecrit  uV/lms  voltage  change  used  as  a  criterion  for 
c  identifying  spikes  in  EEG  channels 

c  voltcrit  Absolute  voltage  for  identifying  an  outlying  data  point 

c  rmscrit  The  critical  epochwise^rms  for  identifying  a  noisy  channel 

c  deadcrit  The  minimum  epochwise  rms  for  identifying  a  dead  EEG  channel 

c  samprate  The  integer  EEG  sampling  rate  (points  /  sec) 

c 

c  The  file  is  in  the  format  of  one  variable  per  line,  free  format, 
c 

c  Important  Program  Variables: 
c 

c  vl(p)  -  Voltage  at  time  point  p,  of  a  vector  of  EEG  data  (or  sometimes 

c  EOG  data).  The  vector  contains  real  values  in  microvolts, 

c  The  points  are  organized  in  temporal  order  from  the  first  point 

c  of  the  epoch  to  the  last  point  of  the  epoch, 

c  v2(p)  -  Voltage  at  point  p  of  the  vertical  EOG.  Organization 
c  of  the  vector  is  like  that  of  vl . 

c  v3(p)  -  Same  as  v2  for  the  lateral  eog. 

c  nptot  -  The  total  number  of  data  points  per  channel  (-  np  *  ne) . 

c  nc  -  #  of  channels, 

c  ne  -  Epochs  of  data  per  channel, 

c  np  -  Points/epoch/channel . 

c  LEYECHAN  -  Integer  channel  number  of  the  lateral  EOG  data, 

c  VEYECHAN  -  Integer  channel  number  of  the  vertical  EOG  data, 

c  nptot  -  The  number  of  epochs  actually  analyzed  *  points/epoch . 
c 

c  Input  Data  File  Structure: 
c 

c  The  EEG  and  EOG  data  files  should  be  direct-access  files  containing 

c  a  vector  of  two-byte  integer  (integer*2)  values  in  each  record.  The 

c  record  organization  is: 
c 

c  epoch  1  channel  1  vector 

c  epoch  1  channel  2  vector 

c 
c 
c 

c  epoch  1  channel  n  vector 

c  epoch  2  channel  1  vector 

c  epoch  2  channel  2  vector 

c 
c 
c 

c  epoch  N  channel  n  vector 

c 

c  Notes: 
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c 

c  The  number  of  channels  cannot  exceed  20,  the  maximum  number  of  epochs 

c  cannot  exceed  700  and  the  maximum  points  per  epoch  cannot  exceed  400. 

c  These  are  limitations  due  to  the  largest  arrays  we  can  run  on  our  computer 

c  and  could  be  increased  on  other  machines  with  more  memory, 

c 

c  Variable  declarations: 
c 

implicit  none 

integer  MAXCHANNELS , MAXEPOCHS 
integer  NETOT , MAXPOINTS ,  VECTORS IZE 
integer  ARTFILE , INFILE , OUTFILE , PARMFILE 
integer  INFOFILE 
integer  STDIN,  STDOUT,  STDERR 


c 


parameter 

(STDIN 

- 

5) 

!  Standard  input 

parameter 

(STDOUT 

- 

6) 

!  Standard  output 

parameter 

(STDERR 

- 

6) 

!  Standard  error  (usually  0) 

parameter 

(MAXCHANNELS 

- 

15) 

!  Maximum  number  of  channels 

parameter 

(MAXEPOCHS 

- 

700) 

!  Maximum  number  of  epochs 

parameter 

(NETOT 

- 

700) 

!  Ditto 

parameter 

(MAXPOINTS 

Mk 

400) 

!  Maximum  number  of  points 

c 

c  May  have  to  use 

some  sort  of  dynamic 

scheme 

if  the  following  get  too  big 

c  (or  just  run  the  1/2  epochs  at  a  time,  which  is  fine  for  an  epochwise 
c  algorithm) 
c 

parameter  (VECTORSIZE 
c 

parameter  (ARTFILE  - 

parameter  (INFILE  - 

parameter  (OUTFILE  - 

parameter  (PARMFILE  - 

parameter  (INFOFILE  - 

c 

c  Switches 
c 

logical  verbose, 
filter, 

eegchan (MAXCHANNELS) 


MAXPOINTS  *  MAXEPOCHS) 

10) 

ID 

12) 

13) 

14) 


!  Verbosity  of  displayed  output 
!  Subject  the  data  to  3  point  filter 
!  Indicates  if  the  channel  is  to  be 
!  artifa' t  filtered 


c 


integer 

i.j ,R. 

nn, 
ios , 
reclen , 
recno , 
VEYECHAN , 
LEYECHAN , 
nctot , 
nc , 
np, 

ichan, 


!  Statistics  counter 
!  General  purpose  error  return 
!  Epochchan  record  length  in  bytes 
!  Current  record  number 
!  Vertical  eye  channel 
!  Lateral  eye  channel 
!  Total  number  of  channels 
!  Ditto,  in  another  context 
!  Number  of  points/epoch 
!  Channel  index 


13 


c 


c 


c 


c 


c 


iepoch, 
ipoint , 
nepochs , 
firstepoch, 
lastepoch, 
delta , 
n, 

ival , 
getcwd, 
retval , 
nominalf irst , 
nominallast , 
nominale , 
nptot , 
pass , 

f irstpoint , 
lastpoint, 
ptemp , 

lencurrentpath , 

leninpath, 

lenparmpa  tli , 

leninfopath, 

lenartpath, 

lenoutpath 


Epoch  index 

Sampled  point  index 

Number  of  epochs  in  a  run 

First  epoch  to  process 

Last  epoch  to  process 

Actual  sampling  increment  used 

General  purpose  counter 

Integer  command  line  argument  value 

System  call  to  get  current  working  directory 

General  purpose  return  value 

Nominal  (relative)  first  epoch  number 

Nominal  (relative)  last  epoch  number 

Nominal  (relative)  current  epoch  number 

Total  number  of  points  in  a  sample 

Used  for  controlling  artifact  filtering  loop 

First  point  to  use  in  a  v_()  array 

Last  point  to  use  in  a  v_()  array 

Used  to  sequentially  index  through 

v_()  arrays 

Length  of  constructed  current  path 
Length  of  constructed  input  file  path 
Length  of  output  parameter  file  path 
Length  of  constructed  information  file  path 
Length  of  artifact  file  path 
Length  of  constructed  output  file  path 


integer 

blink(MAXCHANNELS, MAXEPOCHS) , 
voltart(MAXCHANNELS, MAXEPOCHS) , 
clipart(MAXCHANNELS .MAXEPOCHS) , 
rmsart(MAXCHANNELS .MAXEPOCHS) , 
deadchan (MAXCHANNELS .MAXEPOCHS) , 
spike (MAXCHANNELS , MAXEPOCHS) , 
sumart (MAXCHANNELS .MAXEPOCHS) 


!  Blink  indicator 
!  Absolute  voltage  indicator 
!  Voltage  clipping  indicator 
!  Excessive  RMS  indicator 
!  Not  enough  RMS  indicator 
!  Voltage  spike  indicator 
!  Artifact  summary  variable 


integer*2 

dummy (MAXPOINTS ) , 
arts (10) 


!  Used  to  read  input  data 
!  Artifact  indicator  array 


real 

vl (VECTORSIZE) , 
v2 (VECTORS I ZE) , 
v3 (VECTORSIZE) , 
r dummy (MAXPOINTS) 


!  Vector  used  to  hold  EEG 
!  Lateral  EOG  vector 
1  Vertical  EOG  vector 
!  Input  data  converted  to  float 


real 

rms (MAXCHANNELS, MAXEPOCHS ) ,  !  RMS  (per  channel  and  epoch) 

me an (MAXCHANNELS .MAXEPOCHS)  !  Mean  voltage  (per  channel  and  epoch) 


real 

voltcrit , 
blinkcrit , 
spikecrit , 
rmscrit , 
deadcrit , 


!  Absolute  voltage  criterion 
!  Blink  rate  of  change  criterion 
!  Spike  rate  of  change  criterion 
!  Excessive  RMS  criterion 
!  Dead  channel  criterion 
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clipcrit, 


Clipping  criterion 


.  clipcrit2,  ! 

cov ,  ! 

ss ,  t 

.  tx,  ! 

ty,  ! 

,  vladj ,  ! 

v2adj ,  I 

rval ,  ! 

b(MAXCHANNELS) ,  ! 

.  samprate,  ! 

d,  ! 

.  realpoints ,  ! 

.  realint,  f 

.  realn  ! 

c 

character*256 

arg,  ! 

.  cval  ! 


Adjusted  eog  clipping  criterion 

Calculated  covariance 

Sum  of  sc  ares 

Sum  of  x 

Sum  of  y 

Value  of  vl  with  any  corrections  applied 
Value  of  v2  with  any  corrections  applied 
Real  command  line  value 

Slope  of  correction  line  for  a  given  channel 
Sampling  rate 
Temporary  variable 
Float(nptot) 

Actual  sampling  interval 

Float(number  of  points  checked  per  channel) 


Command  line  string 
Character  command  line  value 


c 


c 


c 

c 

c 


character*l 

option,  ! 

optiontype  ! 

character*256 

currentpath,  ! 

artpath,  ! 

inpath,  ! 

outpath ,  ! 

parmpath ,  1 

infopath  ! 


"transparm"  returned  option  letter 
"transparm"  returned  option  type 


Current  directory  path 
Artifact  file  path  (artifacts) 

Input  file  path  (epochchan) 

Output  file  path  (epochchan. c) 
Parameter  file  path  (artifact. p2) 
Path  to  information  file  (epochinfo) 


Initialize  command  line  variables  with  default  values 


n  -  0 
nc  -  0 
np  -  0 

currentpath  -  '  ' 

inpath  -  'epochchan' 

call  removespaces  (inpath , leninpath) 

outpath  -  ' epochchan. c ' 

call  removespaces  (outpath, lenoutpath) 

parmpath  -  'artifact. p2 ' 

call  removespaces  (parmpath , lenparmpath) 

artpath  -  'artifacts' 

call  removespaces  (artpath , lenartpath) 
VEYECHAN  -  1 
LEYECHAN  -  2 
filter  -  .false, 
do  i-  1,  MAXCHANNELS 
eegchan(i)  -  .true, 
enddo 

firstepoch  -  1 
lastepoch  -  10 
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verbose  -  .false. 


c 

c  Process  command  line  arguments 
c 

1001  n  -  n  +  1 

call  getarg  (n,argN 
if  (arg.eq.'  ')  goto  1002 

call  transparm  (arg,  option,  optiontype  ,  ival ,  rva1. , cval) 
c 

if  ( option. eq. 'c' )  then 
nc  -  ival 
nctot  -  nc 

else  if  (option,  eq . 'v' )  then 
VEYECHAN  -  ival 
else  if  (option. eq .' 1 ' )  then 
LEYECHAN  -  ival 
else  if  (option. eq. 'F' )  then 
firstepoch  -  ival 
else  if  (option. eq. 'L' )  then 
lastepoch  -  ival 
else  if  (option. eq. 'e' )  then 
eegchan(ival)  -  .false, 
else  if  (option. eq . 'm' )  then 
filter  -  .true, 
else  if  (option. eq. 'V' )  then 
verbose  -  . true . 
else  if  (option. eq. 'I')  then 
inpath  -  cval 

call  removespaces  (inpath, leninpath) 
else  if  (option. eq. 'O' )  then 
outpath  -  cval 

call  removespaces  (outpath, lenoutpath) 
else  if  (option. eq. 'P' )  then 
parmpath  -  cval 

call  removespaces  (parmpath, lenparmpath) 
else  if  (option. eq. 'A' )  then 
artpath  »  cval 

call  removespaces  (artpath , lenartpath) 
endif 
goto  1001 

1002  continue 
c 

c  If  no  parameters,  display  help  and  exit 
c 

if  (n.eq.l)  then 

write  (STDOUT ,  *) 
write  (STDOUT,*) 
write  (STDOUT,*) 
write  (STDOUT,*) 
write  (STDOUT,*) 
write  (STDOUT,*) 
write  (STDOUT,*) 
write  (STDOUT,*) 
write  (STDOUT,*) 
write  (STDOUT,*) 


'artfil:  Evoked  potential  artifact  detector' 

'  and  filter  --  Single  epoch  corrections' 

'  -cCnumber  of  channels>,  required' 

'  [-vCvertical  eye  channel>] ,  default  1' 

'  [-lclateral  eye  channei>] ,  default  2' 

'  [-F<first  epoch  to  process>] ,  default  1' 

'  [-Ldast  epoch  to  process>]  ,  default  10' 

'  (-e(channel  to  omit  from  artifact  filtering)],' 
'  default  no  omissions' 


16 


write  (STDOUT,*)  ' 
write  (STDOUT,*)  ' 
write  (STDOUT,*)  ' 
write  (STDOUT,*)  * 
write  (STDOUT,*)  ' 
write  (STDOUT,*)  ' 
write  (STDOUT,*)  ' 
write  (STDOUT,*)  ' 
write  (STDOUT,*)  ' 
write  (STDOUT,*)  ' 
write  (STDOUT,*)  ' 
goto  32767 
endif 


[-m(perform  3  point  filtering  on  data)],' 
default  no  filtering' 

( -V(verbose) ] ,  default  no  verbose' 
[-P<artifact  parameter  file  path  narae>] , ' 
default  “artifact. p2" ' 

[-Kinput  file  path  name>] , ' 
default  "epochchan"' 

[-0<output  file  path  name>] , ' 
default  "epochchan. c" ' 

[-A<artifact  output  file  path  name>] , ' 
default  "artifacts"' 


c 

c  Setup  the  file  path  information 
c 

retval  -  getcwd  (currentpath) 
call  removespaces  (currentpath , lencurrentpath) 
c 

i  -  0 

if  (nctot.le.O)  then 

write  (STDERR,*)  'artfil:  Number  of  channels  must  not  be  ' 
.  'zero  or  less’ 

i  -  1 
endif 

if  (i.ne.O)  then 
goto  32767 
endif 
c 

infopath  -  'eyochinfo' 

call  removespaces  ( infopath , leninfopath) 
c 

c  Get  number  and  length  of  epochs  from  information  file 
c 

open  (unit-INFOFILE,  file-infopath(l : leninfopath) , 

.  status-'old' ,  iostat-ios) 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artful:  Error  opening  epochinfo,  ', 

.  ' iostat  is  ' , ios 

goto  32767 
endif 

rewind  (INFOFILE) 


read  (INFOFILE,*)  nepochs  !  number  of  epochs 

read  (INFOFILE,*)  np  !  number  of  points  per  epoch 

close  (INFOFILE) 
c 

lastepoch  -  nepochs 
if  (verbose)  then 

write  (STDOUT,*)  'Number  of  epochs:  ’, nepochs 
write  (STDOUT,*)  'Number  of  points/epoch:  ' ,np 
endif 
c 

c  Read  the  artifact  parameter  file: 


17 


open  (PARMFILE,  f ile-parmpath(l : lenparmpath) , 

.  status-'old' ,  iostat-ios) 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  opening  parameter  file,  ', 

' iostat  is  ' , ios 

goto  32767 
endif 

rewind  (PARMFILE) 
read  (PARMFILE,*)  blinkcrit 
read  (PARMFILE,*)  spikecrit 
read  (PARMFILE,*)  voltcrit 
read  (PARMFILE,*)  rmscrit 
read  (FARMFILE,*)  deadcrit 
read  (PARMFILE,*)  clipcrit 
read  (PARMFILE,*)  clipcrit2 
read  (PARMFILE,*)  samprate 
close  (PARMFILE) 

if  (verbose)  then 

write  (STDOUT,*)  'Blink  criteria  is  '.blinkcrit,'  uV/10ms' 
write  (STDOUT,*)  'Spike  criteria  is  '.spikecrit,'  uV/ms ' 
write  (STDOUT,*)  'Volt  criteria  is  '.voltcrit,'  uV' 
write  (STDOUT,*)  'RMS  criteria  is  '.rmscrit,'  uV' 
write  (STDOUT,*)  'Dead  criteria  is  '.deadcrit,'  uV' 
write  (STDOUT,*)  ' EEG  clipping  criteria  is  '.clipcrit,'  uV' 
write  (STDOUT,*)  'EOG  clipping  criteria  is  ',clipcrit2,  '  uV' 
write  (STDOUT,*)  'Sampling  rate  is  '.samprate,'  Hz' 
endif 

if  (verbose)  then 

write  (STDOUT,*)  'Vertical  eye  channel  is  ' .VEYECHAN 
write  (STDOUT,*)  'Lateral  eye  channel  is  '.LEYECHAN 
endif 

if  (verbose)  then 

write  (STDOUT,*)  'Gain  of  vertical  eye  channel  is  ' , 

.  clipcrit2/clipcrit , '  %  re  other  eeg  channels' 

endif 

Initialize  data  arrays: 
do  ichan  -  1 ,  nc 

do  iepoch  -  firstepoch,  lastepoch 
blink( ichan, iepoch)  -  0 
voltart( ichan, iepoch)  -  0 
clipart(ichan, iepoch)  -  0 
rms(ichan, iepoch)  -  0. 
rmsart( ichan, iepoch)  -  0 
mean( ichan , iepoch)  -  0, 
deadchan( ichan, iepoch)  -  0 
spike (ichan, iepoch)  -  0 
enddo 
enddo 


DEFINE  SOME  VARIABLES: 


c  Nominalfirst  and  nominallast  are  nominal  epoch  numbers, 
c  Nominal  first  is  always  1. 

c  This  is  the  first  value  in  the  sequence  bounded  by 
c  firstepoch  and  lastepoch.  nominallast  is  the  ordinal  value 
c  of  the  last  epoch  in  the  sequence.  Nptot  is  the  rank  of  the 
c  vl-v3  vectors, 
c 

nominalfirst  -  1 

nominallast  -  lastepoch  -  firstepoch  +  1 
nptot  -  np  *  (lastepoch  -  firstepoch  +  1) 

if  (nptot. gt.VECTORSIZE)  then 

write  (STDERR,*)  'artfil:  Too  many  epochs,  rerun  with  fewer' 
goto  32767 
endif 
c 

reclen  -  np  *  2 
c 

if  (verbose)  then 

write  (STDOUT,*)  'First,  last  epochs:  ', 

.  firstepoch,  lastepoch 

write  (STDOUT,*)  'Nominal  first,  last  epochs:  ', 

.  nominalfirst,  nominallast 

write  (STDOUT,*)  'Number  of  points  total  (nptot):  ', 

,  nptot 

write  (STDOUT,*)  'Number  of  points  allowed:  ', 

VECTORS I ZE 

write  (STDOUT,*)  'Number  of  points/epoch:  ', 
np 

endif 

c 

c  Correct  blink  and  spike  voltage -change  criteria  for 
c  quantization  errors  due  to  the  obtained  intervals 
c  between  time  points, 
c 

realpoints  -  .01  *  float (samp rate) 
c 

c  Set  the  calculating  interval  in  points  to  the  value  nearest 
c  the  number  of  points  giving  10  ms. 
c 

delta  -  realpoints  +0.5 
if  (delta. It. 1)  then 
delta  -  1 
endif 

realint  -  float(delta)/float(samprate)  !  gives  seconds 

!  corresponding  to  delta 

blinkcrit  -  blinkcrit  *  realint  /  .01  !  adjust  blinkcrit  by 

!  ratio  of  delta  to  10  ms 
spikecrit  -  spikecrit  *  1000 . /float(samprate) 
if  (verbose)  then 

write  (STDOUT,*)  'Adjusted  blink  criteria  is  ', 
blinkcrit,'  uV/sample__interval ' 
write  (STDOUT,*)  'Adjusted  spike  criteria  is  ' , 
spikecrit,'  uV/pt' 

endif 
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c 

c  Now,  calculate  "realn",  the  number  of  points  to  be  checked  per  channel, 
c  expressed  as  a  floating  point  number 
c 

realn  -  float (np  -  2  *  delta) 
if  (verbose)  then 

write  (STDOUT,*)  'The  number  of  points  checked  per  chan,  is  ' , 

.  realn 

endif 
c 

c  Open  the  input  and  output  files 
c 

open  (INFILE,  f ile-inpath(l : lenlnpath) ,  status-'old' , 

.  access-'dir',  recl-reclen,  iostat-ios) 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  opening  input  file,  ', 

' iostat  is  ' , ios 

goto  32767 
endif 

open  (OUTFILE,  file-outpath(l : lenoutpath) ,  status- ' unk ' , 

.  access-'dir',  reci-reclen) 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  opening  output  file,  ', 

.  'iostat  is  ',ios 

goto  32767 
endif 
c 

c  We  need  the  eye  movement  data  early  on,  so  we 
c  get  the  vertical  and  horizontal  eye  data  now 

r, 

c  First,  the  vertical  eye  channel... 
c 

if  (verbose)  then 

write  (STDOUT,*)  'Reading  vertical  eye  channel  data. . . ' 
endif 

ichan  -  VEYECHAN 
ptemp  -  0 

do  iepoch  -  firstepo^h,  lastepoch 

recno  -  nctot  *  (iepoch  -  1)  +  ichan 

call  quickio  (INFILE,  inpath(l : leninpath) ,  'old',  'read',  recno, 

.  dummy,  np,  ios) 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  reading  input  file,  ', 

.  ' iostat  is  ' , ios 

goto  32767 
endif 
c 

c  Note  that  the  filtering  subprogram  "smooth3p"  scales  everything  down  by  a 
c  factor  of  10,  to  account  for  the  assumption  that  the  input  file  is  in  1/10  uV 
c.  units  (i.e.  10.6  uV  is  represented  by  106  in  the  file).  If  filtering  is  not 
c  performed,  then  we  need  to  adjust  the  values  ourselves, 
c 

if  (filter)  then 

call  smoo th3p( dummy, rdummy.np) 

else 
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do  1  -  1,  np 

rdummy(i)  -  dummy(i)  /  10.0 
enddo 
end  if 

do  ipoint  -  1,  np 

ptemp  -  ptemp  +  1 
v2(ptemp)  -  rdummy( ipoint) 
et.ddo 
enddo 
c 

c  And  then  the  lateral  eye  channel 
c 

if  (verbose)  then 

write  (STDOUT,*)  'Reading  lateral  eye  channel  data...' 
endif 

ichan  -  LEYECHAN 
ptemp  -  0 

do  iepoch  -  firstepoch,  lastepoch 

recno  -  nctot  *  (iepoch  -  1)  +  ichan 

call  quickio  (INFILE,  inpath(l : leninpath) ,  'old',  'read', 

.  recno,  dummy,  np,  ios) 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  reading  lateral  eye  ', 

.  'channel  data,  iostat  is  ' , ios 

goto  32767 
endif 

if  (filter)  then 

call  smooth3p ( dummy , r dummy ,np) 

else 

do  i  -  1,  np 

rdummy(i)  -  dummy(i)  /  10.0 
enddo 
endif 

do  ipoint  -  1,  np 

ptemp  -  ptemp  +  1 
v3(ptemp)  -  rdummy( ipoint) 
enddo 
enddo 
c 

c  DO  A  PRELIMINARY  CHECK  OF  THE  VERTICAL  EOG : 
c 

c  Before  doing  anything  else,  check  the  vertical  eog  channel  for 
c  blinks  and  clipping.  Do  not  check  it  for  spikes,  in  order  not  to 
c  confuse  blinks  with  spikes, 
c 

if  (verbose)  then 

write  (STDOUT,*)  'Checking  vertical  EOG  for  blinks  and  clipping. . 
endif 

iepoch  -  firstepoch  -  1 
do  nominale  -  nominalf irst ,  nominallast 
iepoch  -  iepoch  +  1 

j  -  o 

k  -  0 
d  -  0.0 

firstpoint  -  (nominale  -  1)  *  np  +  1 
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lastpoint  -  firstpotnt  +  np  -  1 
do  ipoint  -  firstpoint,  lastpoint-delta 

if  (abs(v2(ipoint) -v2(ipoint+delta)) . gt.  blinkcrit)  then 
if  (j.eq.O)  then 
j  -  ipoint 

<3  -  abs(v2(ipoint) -v2(  ipoint+delta)  ) 
endif 

blink(VEYECHAN, iepoch)  -  1 
endif 

if  (abs(v2(ipoint)) . gt . clipcrit2)  then 
if  (k.eq.O)  then 
k  -  ipoint 
endif  , 

clipart(VEYECHAN , iepoch)  -  1 
endif 
enddo 

if  (verbose . and. blink(VEYECHAN , iepoch) .eq . 1)  then 
print  *, 'Epoch  '.nominale,'  blink  starting  at  ' , 

.  j -firstpoint , '  Vdiff-',d 

endif 

if  (verbose .and. clipart(VEYECHAN, iepoch) .eq. 1)  then 
print  *, 'Epoch  '.nominale,'  clipped  at  pt  ', 

.  k-firstpoint,  '  voltage  is  ',v2(k) 

endif 
enddo 
c 

c  First,  write  vertical  EOG  directly  to  output  file 
c 

if  (verbose)  then 

write  (STDOUT ,*)  'Writing  Veog  to  output...' 
endif 

ichan  -  VEYECHAN 

do  iepoch  -  firstepoch,  lastepoch 

recno  -  nctot  *  (iepoch  -  1)  +  ichan 

call  quickio  (INFILE,  inpath(l : leninpath) ,  'old',  'read',  recno, 
dummy,  np,  ios) 
if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  reading  vertical  EOG  ', 

.  'data,  iostat  is  ' , ios 

goto  32767 
endif 

call  quickio  (OUTFILE,  outpath( 1 : lenoutpath) , 'unk' , 'write ', recno , 
dummy,  np,  ios) 
if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  writing  vertical  EOc  ', 

'data,  iostat  is  ' , ios 

goto  32767 
endif 
enddo 
c 

c  Now,  process  remaining  channels  (Lateral  eog  and  all  eeg) 
c 

if  (verbose)  then 

write  (STDOUT,*)  'Processing  remaining  channels . . . ' 
endif 
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c 

pass  -  0 
ichan  -  0 
c 

150  continue.  !  Top  of  loop  through  channels 

c 

c  First,  check  and  correct  the  lateral  eye  channel.  From  that  point  on, 
c  process  the  eeg  channels  in  order.  The  vertical  eye  channel  is  not 
c  checked, 
c 

if  (pass.eq.O)  then 
ichan  -  LEYECHAN 
pass  -  pass  +  1 
else  if  (pass.eq.l)  then 
ichan  -  3 
pass  -  pass  +  1 

else 

ichan  -  ichan  +  1 
pass  -  pass  +  1 
endif 
c 

if  ( ichan. gt .nc)  then 
goto  3100 
endif 
c 

c  Skip  the  non-eeg  channels,  but  write  the  data  to  the  output  file 
c  unchanged 
c 

if  ( .not. eegc.han( ichan) )  then 

do  iepoch  -  firs iepoch,  lastepoch 

recno  -  nctot  *  (iepoch  -  1)  +  ichan 

call  quickio  (INFILE,  inpath(l : leninpath) ,  'old',  'read', 

.  recno,  dummy,  np,  ios) 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  reading  non-EEG  ', 
'data,  iostat  is  ' , ios 

goto  32767 
endif 

call  quickio  (OUTFILE,  outpath(l : lenoutpath) , 'unk ', 'write ' . 
.  recno,  dummy,  np,  ios) 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  writing  non-EEG  ', 
'data,  iostat  is  ' , ios 

goto  32767 
endif 
enddo 
goto  150 
endif 


c 

c  If  we  are  analyzing  EOG  data,  we  read  the  EOG  into 
c  vl  from  v3 .  Otherwise  we  get  the  EEG  from  file, 
c  Remember:  vl  is  the  work  vector,  v2  is  the  Veog  vector, 
■z  and  v3  is  the  Leog  vector 


c 

if  ( ichan. eq.VEYECHAN)  then 
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do  i  -  1,  nptot 
vl(i)  -  v2 (i) 
enddo 

else  if  ( ichan . eq . LEYECHAN)  then 
do  1  -  1,  nptot 
vl(i)  -  v3(i) 
enddo 

else 

ptemp  -  0 

do  iepoch  -  firstepoch,  lastepoch 

recno  -  nctot  *  (iepoch  -  1)  +  ichan 

call  quickio  (INFILE,  inpath(l : leninpath) ,  'old', 

'read',  recno,  dunjmy,  np ,  ios) 
if  (ios.ne.O)  then 

write  (STDERR , *)  'artfil:  Error  reading  EEG  ', 

'data,  iostat  is  ' , ios 

goto  32767 
endif 

if  (filter)  then 

call  smooth3p ( dummy ,rdummy,np) 

else 

do  i  -  1 ,  np 

rdummy(i)  -  dummy (i)  /  10.0 
enddo 
endif 

do  ipoint  -  1 ,  np 

ptemp  -  ptemp  +  1 
vl (ptemp)  -  rdummy( ipoint) 

enddo 

enddo 

endif 

c ********************  ********************************************************** 
c 

c  BEGIN  THE  PRE- CORRECTION  ARTIFACT  CHECK, 
c 

c  A  final  check  for  absolute  voltage  artifacts  is  performed  after  lateral 
c  eye  movement  artifacts  are  removed  from  the  eog,  in  the  eye-artifact 
c  filter  section  below 
c 

iepoch  -  firstepoch  -  1 
do  nominale  -  nominalf irst ,  nominallast 
iepoch  -  iepoch  +  1 

firstpoint  -  (nominale  -  1)  *  np  +  1 

lastpoint  -  firstpoint  +  np  -  1 

do  ipoint  -  firstpoint+  delta,  lastpoint-delta 

rms ( ichan , iepoch)  -  rms ( ichan , iepoch)+vl ( ipoi nt)**2 . 
mean( ichan , iepoch)  -  mean( ichan , iepoch)+vl ( ipoint) 
c 

c  Spike  detection: 
c 

if  ( ichan. ne.VEYECHAN)  then 

if  ( (vl( ipoint) -vl (ipoint+1) )  . gt .  spikecrit)  then 
spike( ichan, iepoch)  -  1 
endif 
c 
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c  Voltage  clipping  detection:  Notice  that  there  are  separate  criteria 
c  for  the  eye  and  frontal  channels, 

c 

if  ( ( ichan. eq . VEYECHAN) . or . (ichan. eq . LEYECHAN) )  then 
if  (abs (vl( ipoint) ) . gt . clipcrit2)  then 
clipart(ichan, iepoch)  -  1 
end  if 

else  !  other  channels  (normal  EEG  channels) 

if  (abs(vl(ipoint) ) . gt . clipcrit)  then 
clipart(ichan, iepoch)  -  1 
if  (verbose)  then 

print  *, 'Epoch  '.iepoch,'  channel  '.ichan, 

'  clipping  at  pt  '.ipoint,'  v-' , 

.  vl(ipoint) 

endif 
end  if 
endif 
endif 
enddo 
c 

c  Finish  calculating  epoch  means  and  mean  squared  values: 
c 

rms( ichan, iepoch)  -  ( rms (ichan, iepoch) - 
.  mean( ichan, iepoch) **2 .0  /  realn)  /  realn 

raean( ichan, iepoch)  -  mean( ichan, iepoch)/realn 
c 

c  Try  to  detect  possible  rounding  errors  or  precision  problems.  Do  so 
c  without  stopping  the  program,  but  display  an  error  message  so  the  problem 
c  can  be  investigated, 
c 

if  ( rms (ichan, iepoch) .gt. 0.0)  then 

rms (ichan, iepoch)  -  sqrt(rms(ichan, iepoch)) 

else 

write  (STDERR,*)  'artfil:  RMS  for  epoch  '.iepoch, 

.  'channel  '.ichan, 

'  was  bad  (sumsq  <  sum**2/n) .  ' , 

.  '  Setting  to  0.0. ' 

rms (ichan, iepoch)  -  0. 
endif 
c 

c  Excessive  RMS  detection: 
c 

if  (rms(ichan, iepoch) .gt. rmscrit)  then 
rmsart( ichan, iepoch)  -  1 
endif 
enddo 
c 

c  Dead  channel  detection:  Note  that  eye  channels  are  never  declared  dead, 
c 

if  (( ichan. ne . LEYECHAN) . and .( ichan , ne .VEYECHAN) )  then 
do  iepoch  -  firstepoch,  lastepoch 

if  (rms(ichan, iepoch) .It. deadcrit)  then 
deadchan( ichan, iepoch)  -  1 

endif 

enddo 
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endif 


c 

c  END  PRE- CORRECTION  ARTIFACT  CHECK 
c 

C********************************************'**  *******  *************************** 

c 

c  BEGIN  EYE  MOVEMENT  FILTER  AND  FINAL  ARTIFACT  CHECK, 
c 

c  The  eeg  is  checked  for  absolute  amplitude  artifacts  at  the  end  of  this 
c  block  of  code,  after  subtracting  the  scaled  lateral  eog  from  the  eeg 
c 

c  Blink  correction: 

c  j 

c  Regress  the  eeg  on  the  eog  in  regions  where  the  eog  is  changing  rapidly 
c  and  get  the  proportionality  constants,  b(c),  over  all  channels.  Epochs 
c  with  spikes  in  the  eeg  or  clipping  in  the  eog  or  eeg  are  not  used, 
c 

c  An  inefficiency  of  this  section  of  code  is  that  blinks  are  detected  each 
c  time  the  code  is  executed, 
c 

c  The  lateral  eog  is  corrected  first, 
c 

if  (ichan.ne .VEYECHAN)  then 
if  (verbose)  then 

write  (STDOUT , *)  »  Regressing  EEG  on  EOG...' 
endif 
c 

c  Process  each  epoch  separately 
c 

iepoch  -  firstepoch  -  1 
do  nominale  -  nominalf irst ,  nominallast 
iepoch  -  iepoch  +  1 

if  (  (blink(VEYECHAN, iepoch)  .eq.  1)  .and. 

( spike ( ichan, iepoch)  .eq.  0)  .and. 

.  (clipart(VEYECHAN , iepoch)  .eq.  0)  .and. 

(clipart(ichan, iepoch)  .eq.  0))  then 

c 

c  Place  the  pointers  at  the  proper  place  in  the  vector 
c 

firstpoint  -  (nominale  -  1)  *  np  +  1 
lastpoint  -  (firstpoint  +  np)  -  1 
c 

c  Form  the  regression  sums  for  this  epoch 
c 

cov  -  0. 
ss  -  0 . 
tx  -  0. 
ty  -  0. 
nn  -  0 

do  ipoint  -  f irstpoint+delta ,  lastpoint- delta 

if  (abs(v2(ipoint+delta) -v2 ( ipoint -delta) ) . gt . 
blinkcrit)  then 
c 

vladj  -  vl(ipoint) 
v2adj  -  v2( ipoint) 
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tx  -  tx  +  v2adj 
ty  -  ty  +  vladj 
nr.  -  nn  +  1 

cov  -  cov  +  vladj  *  v2adj 
ss  -  ss  +  v2adj  *  v2adj 
end  if 
enddo 


c 

c  Calculate  the  blink  scaling  constants  for  each  channel 
c 

b( ichan)  -  0.0 
if  (ss.eq.O.)  then 

if  (nn.gt.O)  then 
write  (STDERR,*) 

'artfil:  Attempt  to  divide  VRTEYE 
.  'covariance  by  zero  at  channel  ' , ichan 

endif 

else 

if  (nn.gt.O)  then 

ss  -  ss  -  (tx**2 . )/float(nn) 
cov  -  cov  -  (tx*ty)/float(nn) 
if  (ss.ne.0.0)  then 

b(ichan)  -  cov  /  ss 
endif 
endif 
endif 
c 

c  Apply  the  correction  factor.  If  something  went  wrong  with  the  calculation 
c  of  b(),  then  It  was  set  to  zero,  causing  no  change  to  the  data, 
c 

do  ipoint  -  firstpoint,  lastpoint 

vl(ipcint)  -  vl(ipoint)  -  b(ichan)*v2( ipoint) 
enddo 

endif  !  blink 

enddo  !  nominale 

endif  !  .not.  VEYECHAN 

c 

c  End  blink  correction 
c 

c  If  the  Leog  was  just  blink  corrected,  replace  the  uncorrected  data 
c  in  V3  with  corrected  data, 
c 

if  (ichan  .eq.  LEYECHAN)  then 
do  i  -  1 ,  nptot 
v3(i)  -  vl(i) 
enddo 
endif 
c 

c  Lateral  eye  movement  correction: 
c 

c  The  lateral  eye  channel  is  corrected  for  blinks  on  the  first  pass.  The 
c  blink  corrected  lateral  eog  channel  is  used  to  correct  the  eeg 
c  thereafter.  Once  again,  epochs  with  spikes  or  clipping  are  not  used, 
c  Also,  once  again,  the  lateral  eog  is  not  corrected  for  lateral  eye  movement. 
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c 


c 


if  ( ichan . ne . LEYECHAN)  then 
cov  -  0. 
ss  -  0 . 
tx  -  0 . 
ty  -  0. 
nn  -  0 

iepoch  -  firstepoch  -  1 
do  nominale  -  nominalf irst ,  nominallast 
iepoch  -  iepoch  +  1 

if  ( (spike (LEYECHAN, iepoch)  .eq.  0)  .and. 

(spike( ichan, iepoch)  .eq.  0)  .and. 

(clipart (LEYECHAN /iepoch)  .eq.  0)  .and. 

(clipart(ichan, iepoch)  .eq.  0))  then 

firstpoint  -  (nominale  -  1)  *  np  +  1 
lastpoint  -  firstpoint  +  np  -  1 
do  ipoint-f irstpoint+delta ,  lastpoint-delta 
tx  -  tx  +  v3(ipoint) 
ty  -  ty  +  vl(ipoint) 
nn  -  nn  +  1 

cov  -  cov  +  vl(ipoint)  *  v3(ipoint) 
ss  -  ss  +  v3(ipoint)  *  v3(ipoint) 
enddo 
endif 
enddo 

if  (ss.eq.0.)  then 

if  (nn.gt.0)  then 

write  (STDERR , *)  'artfil:  Attempt  to  divide  LATEYE 
'covariance  by  zero  at  channel  ' ,  ichan 

endif 

b( ichan)  -  0.0 

else 

if  (nn.gt.0)  then 

ss  -  ss  -  ( tx**2 . )/float(nn) 
cov  -  cov  -  (tx*ty)/float(nn) 
b(ichan)  -  cov  /  ss 

else 

b( ichan)  -  0. 
endif 
endif 


c 

c  Subtract,  the  scaled  lateral  eog  and  check  the  resulting  data 
c  for  absolute  amplitude  artifacts.  Epochs  with  spikes  or  clipping 
c  are  not  corrected.  Also,  the  "delta"  points  at  the  beginning 
c  and  end  of  each  vector  are  not  checked  for  artifacts.  The  lateral 
c  eog  is  not  checked  for  amplitude  artifacts, 
c 

if  (abs(b(ichan) )  , gt .  0.)  then 
iepoch  -  firstepoch  -  1 
do  nominale  -  nominalf irst ,  nominallast 
iepoch  -  iepoch  +  1 
c 

c  Check  for  identified  artifacts: 
c 
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if  ( (spike (LEYECHAN, iepoch)  .eq.  0)  .and. 

( spike ( ichan, iepoch)  .eq.  0)  .and. 

(clipart (LEYECHAN, iepoch)  .eq.  0)  .and. 

(clipart(ichan, iepoch)  .eq.  0))  then 

c 

c  Subtract  the  scaled  eog  if  the  epoch  is  good: 
c 

firstpoint  -  (nominale-l)*np+l 
lastpoint  -  firstpoint+np-1 
do  ipoint  -  firstpoint,  lastpoint 

vl(ipoint)  -  vl(ipoint) -b(ichan)*v3(ipoint) 


c  Check  the  corrected  eeg  for  large  absolute  voltages: 
c 

if  (ichan  ,ne.  LEYECHAN)  then 

do  ipoint  -  f irstpoint+delta ,  lastpoint-delta 
if  (abs(vl(ipoint)) . gt . voltcrit)  then 
voltart(ichan, iepoch)  -  1 
endif 
enddo 

endif  !  .not.  lateral  eye  channel 

endif  !  no  spike  artifacts 

enddo  !  epoch 

endif 
endif 


c 

c  END  LATERAL  EOG  CORRECTION 


c 

c  END  EYE  MOVEMENT  FILTER 


c 

c  Convert  to  integer*2,  rescale  (converts  to  uV) ,  and  write 
c  corrected  EEG  to  the  output  file 
c 

if  (verbose)  then 

print  *, 'Writing  all  epochs  for  channel  '.ichan 
endif 
ptemp  -  0 

do  iepoch  -  firstepoch,  lastepoch 
do  ipoint  -  1 ,  np 

ptemp  -  ptemp  +  1 
c 

c  The  following  results  are  rounded  before  truncating 
c 

if  (vl(ptemp)  . ge .  0)  then 

dummy (ipoint)  -  vl(ptemp)  *  10.  +  .5 

else 

dummy ( ipoint)  -  vl(ptemp)  *  10.  -  .5 
endif 
enddo 


recno  -  nctot  *  (iepoch- 1)  +  ichan 

call  quickio  (OUTFILE,  outpath( 1 : lenoutpath) ,  'unk', 
'write',  recno,  dummy,  np,  ios) 
if  (ios.ne.0)  then 
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write  (STDERR,*)  'artfil:  Error  writing  to  output  ', 

.  'file,  iostat  is  ' , ios 

goto  32767 
end  if 
enddo 
c 

c  End  of  channel  loop 
c 

goto  150  !  process  next  channel 

c 

c  Finished  with  channels,  clean  up 
c 

3100  continue 
c 

c  Fill  out  the  final  records  for  the  last  epoch 
c 

do  ipoint  -  1,  np 

dummy (ipoint)  -  0 
enddo 
c 

recno  -  nctot  *  lastepoch 

call  quickio  (OUTFILE,  outpath( 1 : lenoutpath) ,  'unk',  'write',  recno, 

,  dummy,  np ,  ios) 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Err"*  writing  cleanup  records,  ', 

.  'iostat  is  ' , ios 

goto  32767 
endif 
c 

c  Close  the  data  files 

G 

close  (INFILE) 
close  (OUTFILE) 
c 

c  Consolidate  some  of  the  artifacts  into  a  summary  variable  to  be  used 
c  for  deleting  epochs  in  subsequent  processing.  Channels  within  epochs 
c  are  deletable  if  they  contain  clipping,  large  absolute  or  rms  voltages 
c  in  any  channel  except  the  vertical  eog  channel.  If  there  is  clipping 
c  in  the  vertical  or  lateral  eog,  the  entire  epoch  is  marked  for  deletion, 
c 

if  (verbose)  then 

print  *, 'Creating  artifact  summary  file' 
endif 

do  ichan  -  1 ,  nc 

do  iepoch  -  firstepoch,  lastepoch 

if  (( ichan. eq. LEYECHAN) . and. (clipart( ichan, iepoch) . eq . 1) )  then 
if  (verbose)  then 

print  *, 'Epoch  ', iepoch,'  Chan,  '.ichan, 

.  '  -  Lateral  eye  channel  clipping  detected' 

endif 

sumart( ichan, iepoch)  -  1 
endif 
c 

if  ( (ichan. eq.VEYECHAN) . and. (clipart( ichan, iepoch) . eq . 1) )  then 
if  (verbose)  then 
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print  *, 'Epoch  '.iepoch,'  Chan,  '.ichan, 

'  -  Vertical  eye  channel  clipping  detected' 

endif 

sumart( ichan, iepoch)  -  1 
endif 

if  ( (ichan. eq . LEYECHAN) . and. (spike (ichan, iepoch) . eq . 1) )  then 
if  (verbose)  then 

print  *, 'Epoch  '.iepoch,'  Chan,  '.ichan, 

'  -  Lateral  eye  channel  spike  detected' 
endif 

3umart( ichan, iepoch)  -  1 
endif 

if  ( ( ichan. ne . LEYECHAN) . and. ( ichan. ne . VEYECHAN) ) then 
if  (spike(ichan, iepoch) .eq. 1)  then 
if  (verbose)  then 

print  *, 'Epoch  '.iepoch,'  Chan,  '.ichan, 

'  -  Spike  detected' 

endif 

suraart( ichan, iepoch)  -  1 
endif 

if  (voltart(ichan, iepoch) . eq . 1)  then 
if  (verbose)  then 

print  *, 'Epoch  '.iepoch,'  Chan,  '.ichan, 

'  -  Voltage  artifact  detected' 

endif 

sumart( ichan, iepoch)  -  1 
endif 

if  (clipart(ichan, iepoch) . eq . 1)  then 
if  (verbose)  then 

print  *, 'Epoch  '.iepoch,'  Chan,  '.ichan, 

'  -  Clipping  detected' 

endif 

sumart( ichan, iepoch)  «  1 
endif 

if  (rmsart(ichan, iepoch) .eq. 1)  then 
if  (verbose)  then 

print  *, 'Epoch  '.iepoch,'  Chan,  '.ichan, 

'  -  Excessive  RMS  detected' 

endif 

suinart( ichan,  iepoch)  -  1 
endif 

if  (clipart(LEYECHAN , iepoch) . eq . 1)  then 
if  (verbose)  then 

print  *, 'Epoch  '.iepoch,'  Chan,  '.ichan, 

'  -  Artifact  due  to  lateral  eye  clipping' 

endif 

sumart( ichan, iepoch)  -  1 
endif 

if  (clipart (VEYECHAN , iepoch) . eq . 1) then 
if  (verbose)  then 

print  *, 'Epoch  '.iepoch,'  Chan,  '.ichan, 

'  -  Artifact  due  to  vertical  ' , 

'eye  clipping' 
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endif 

sumart( ichan , iepoch)  -  1 
endif 
endif 
enddo 
enddo 
c 

c  Create  the  artifact  summary  file 
c 

artpath  -  'artifacts' 

call  removespaces  (artpath , lenartpath) 
c 

open  (ARTFILE,  f ile-artpath( 1 : lenartpath) ,  status- ' unk' , 

.  access-'dir' ,  recl-20,  iostat-ios) 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  opening  artifact  file,  ', 

.  ' iostat  is  ' , ios 

goto  32767 
endif 
c 

do  iepoch  -  firstepoch,  lastepoch 
do  ichan  -  1 ,  nc 

arts(l)  -  sumart( ichan, iepoch) 

arts(2)  -  iepoch 

arts (3)  -  ichan 

arts(4)  -  blink(ichan, iepoch) 

arts(5)  -  voltart(ichan, iepoch) 

arts(6)  -  rmsart( ichan, iepoch) 

arts(7)  -  deadchan( ichan, iepoch) 

arts(8)  -  spike(ichan, iepoch) 

arts (9)  -  clipart(ichan, iepoch) 

arts (10)  -  rms( ichan, iepoch)  +  .5 

recno  -  (iepoch  -  1)  *  nctot  +  ichan 

write  (ARTFILE,  rec-recno,  iostat-ios)  arts 

if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  writing  to  ', 

.  '  artifact  summary  file,  ', 

'  iostat  is  ' , ios 

goto  32767 
endif 
c 

c  Fill  out  the  records  for  the  last  epoch,  if  neccessary. 
c 

if  (  (iepoch. eq. lastepoch) .and. 

( ichan . eq . nc) . and . (nc . ne , nctot) )  then 
do  i  -  1,  10 

arts(i)  -  0 
enddo 

recno  -  nctot  *  lastepoch 

write  (ARTFILE,  rec-recno,  iostat-ios)  arts 
if  (ios.ne.O)  then 

write  (STDERR,*)  'artfil:  Error  writing  ', 

'cleanup  records  to  artifact  ', 
'summary  file,  iostat  is  ',ios 

goto  32767 
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endif 

endif 

enddo 

enddo 

close  (ARTFILE) 

32767  stop 
end 

************************************************************************** 

subroutine  smooth3p  (dummy,  rdummy,  np) 
c 

c  Converts  integer*2  vectors 
c  and  computes  a  three -point 
c 

integer*2  dummy(np) 
real  rdummy (np) 
integer  np,  p 
real  zl,  z2,  z3,  kl, 
c 

c  The  smoothing  algorithm  is 
c 

kl  -  .025 
k2  -  .05 

z2  -  float (dummy ( 1) ) 
z3  -  float (dummy (2) ) 
c 

do  p  -  2,  np  -  1 
zl  -  z2 
z2  -  z3 
z3  -  float(dummy(p+l) ) 

rdummy(p)  -  kl  *  zl  +  k2  *  z2  +  kl  *  z3 
enddo 
c 

c  Set  the  first  and  last  elements  of  the  smoothed  vector 
c  the  values  of  their  immediate  neighbors, 
c 

rdummy (1)  -  rdummy (2) 
rdummy (np)  -  r dummy (np-1) 
c 

return 

end 

********************************** **************************************** 

subroutine  transparm  (parmstr , retoption, rettype , 
retival , retrval , retcval) 
integer  nnumind, nintind 
parameter  (nnumind  -  15) 

parameter  (nintind  -  12) 

c 

character*(*)  parmstr 
character*256  retcval.parm 
character*l  rettype , retoption , 

numind( nnumind) , intind(nintind) 


to  real,  divides  each  element  by  10, 
smooth . 


k2 


x(t)  -  . 25x( t-2)+ . 5x(t)+ . 25x(t+l) 
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c 


integer  retival , entirelength ,parmlength 
integer  i,j,ios 
real  retrval 


c 


data  numind  /'O', 

I 

data  intind  /'O' , 


c 

c  Given  a  parameter  string  in  the  format  '-'  option  '...string...' 
c 

c  a)  extract  the  option  (one  character) . 

c  b)  determine  the  type  of  the  string  as  integer,  real,  or  character, 

c  c)  convert  integer  and  real  values  to  internal  representation, 

c  d)  return  values  to  calling  program, 

c 

c  Author:  Sam  J.  La  Cour ,  Jr. 
c  Date:  February  20,  1987 

c 

c  First,  get  the  length  of  the  entire  string.  In  order  for  it  to  be  a  valid 
c  option  string,  the  length  must  be  at  least  2. 
c 

call  getlen  (parmstr.i) 
if  (i.lt.2)  then 

retoption  -  '  ' 
rettype  -  's' 
return 
endif 


c 

c  Next,  see  if  it  is  a  valid  option  parameter  string.  In  order  for  this  to 
c  be  true,  the  first  character  must  be  a  followed  by  an  option  character, 

c  followed  by  an  optional  string, 
c 

if  (parmstr(l : 1) ,ne . ' - ' )  then 
retoption  -  '  ' 
rettype  -  'u' 
return 
endif 


c 

c  Ok,  we  have  a  let's  get  the  option  character.  It  can  be  anything 

c  except  a  blank, 
c 

retoption  -  parmstr(2:2) 
c 

c  Now  that  we  have  a  legal  option  character,  the  remainder  of  the  option 
c  string  needs  to  be  parsed,  if  present.  If  the  length  of  the  option 
c  string  is  equal  2,  just  go  ahead  and  exit, 
c 

if  (entirelength , eq . 2)  then 
rettype  -  'n' 
return 
endif 


c 

c  Get  the  actual  length  of  the  rest  of  the  parameter  string, 
c 
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parm  -  parmstr(3:) 
call  getlen  (parm, parmlength) 
if  (parmlength. le .0)  then 
rettype  -  'n' 
return 
endif 
c 

c  Scan  the  parameter  string  for  non-numeric  characters. 

c  For  our  purposes,  numeric  characters  are  'e' , 'E' , ' - ' , and 

c  All  other  characters  classify  the  string  as  non-numeric,  and  therefore 
c  the  type  is  forced  to  be  character.  After  this  is  done,  if  the  the  type  is 
c  determined  to  be  possibly  numeric,  a  further  scan  is  done  to  eliminate  the 
c  'e','E'  and  This  filter  will  determine  if  the  number  is  possibly  an 

c  integer  value.  Once  this  is  done,  internal  reads  will  be  used  to  do  the 
c  actual  conversion  from  string  to  Internal  format.  Any  failures  in  the 
c  conversion,  due  to  things  like  '1.0eel.3-',  will  be  flagged  here.  If  all 
c  is  successful,  the  type  and  value  will  have  been  determined  and 
c  returned  to  the  main  program, 
c 

c  Pass  through  any  numerics,  send  all  others  to  type  CHARACTER 
c 

do  1010  i— 1 , parmlength 
do  1020  j— l,nnumind 

if  (parm(i : i) . eq .numind( j ) )  goto  1010 
1020  continue 

c 

c  A  match  with  the  numeric  values  was  not  found,  the  string  must 
c  be  of  type  CHARACTER,  therefore,  proceed  to  that  section, 
c 

goto  4000 
1010  continue 

c 

c  All  characters  in  the  string  passed  the  numeric  test,  now  lets 
c  try  to  find  out  if  the  number  is  integer  or  real  by  looking  for 
c  characters  which  would  indicate  that  the  value  was  a  real  number, 
c  such  as  decimal  points  or  exponentials, 
c 

c  Pass  through  any  integers,  send  all  others  to  type  REAL 
c 

do  2010  i-1 .parmlength 
do  2020  j-l,nintind 

if  (parm( i : i) . eq . intind( j ) )  goto  2010 
2020  continue 

c 

c  A  match  with  the  integer  values  was  not  found,  the  string  must 
c  be  of  type  REAL,  therefore  proceed  to  that  section, 
c 

goto  3000 
2010  continue 

c 

c  The  value  has  only  legal  integer  characters  in  it,  namely  '0'..'9'  and 
c  Try  to  convert  it.  Any  errors  are  due  to  incorrect  format, 
c 

read  (parm(l : parmlength) ,*, ios tat-ios)  retival 
if  (ios.ne.0)  then 
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rettype  -  'u' 

else 

rettype  -  'i' 
end  if 
return 
c 

c  The  value  has  only  legal  real  characters  in  it,  namely  ' 0 ' . . ' 9 ' , ' e ' , ' E ' , '+ ' , 
c  Try  to  convert  it. 

c 

3000  continue 

read  (parm(l:parmlength) iostat-ios)  retrval 
if  (ios.ne.O)  then 
rettype  -  'u' 

else 

rettype  -  ' r' 
endif 
return 
c 

c  Character  is  the  only  type  left  (a  trash  disposal...) 
c 

4000  continue 

rettype  -  'c' 
retcval  -  parm 
return 
end 

***********iV******************'********'>Sr***'iV,*vk  *********************  ■*•■***■***** 

subroutine  reraovespaces  (instring  newlen) 
c 

c  Removes  spaces  and  nulls  from  a  string  and  determines  the  length  of  the  string 
c 

character*(*)  instring 
character*256  temp 
integer  i,  newlen 
c 

newlen  -  0 
temp  -  '  ' 

do  1000  i-1 , len( instring) 

if  (instring( i ; i) . ne . '  ' . and. instring(i : i) .ne . char(0) )  then 
newlen  -  newlen  +  1 
temp(newlen: newlen)  -  instring( i : i ) 
endif 

1000  continue 

instring  -  temp 

return 

end 

******************************************************************************* 

subroutine  getlen  (str , strlength) 
c 

c  Given  any  length  string,  returns  the  length  of  the  string  with  all 
c  trailing  spaces  removed.  In  other  words,  the  length  of  the  string 
c  scanned  from  the  right  until  the  first  non  blank  character. 
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c 

c  Examples: 
c 

c  'foo  bar'  returns  a  length  of  7 

c  'foo  bar  '  also  returns  a  length  of  7 

c  '  foo  bar'  returns  a  length  of  9 

c 

c  Author :  Sam  J .  La  Cour ,  Jr. 

c  Date:  February  20,  1987 

c 

character*(*)  str 
integer  strlength, i , j 
c 

c  Ignore  the  trivial  case  and  exit 
c 

j  -  len(str) 
if  (j.le.O)  then 
strlength  -  0 
return 
end  if 
c 

c  Scan  from  the  right,  constantly  updating  'strlength'.  Exit  when  the 
c  entire  string  has  been  scanned.  If  no  non -blanks  are  found,  the  length 
c  will  be  zero, 
c 

strlength  -  j 
do  1000  i-1 , j 

if  (str(strlength: strlength) ,ne . '  ')  go  to  32767 
strlength  -  strlength  -  1 
1000  continue 
32767  return 
end 

**********************************************  ***•**************************•*•**-*•■* 

subroutine  casefold  (string) 

c 

c  Routine  to  change  all  characters  in  a  string  to  upper  case  for  comparison 
c  This  is  most  useful  on  things  like  "y"  or  "Y"  inputs,  where  it  doesn't  matter 
c  what  the  case  of  the  response  is. 
c 

character*(*)  string 
integer  i , j ,1 
c 

1  -  len(string) 
do  1000  i-1,1 
j  -  ichar (string( i : i) ) 
if  (j  ge . 97 . and. j . le . 122)  then 
j  “  j  -  ^ 

string(i:i)  -  char(j) 
endif 

1000  continue 
return 
end 
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subroutine  quickio  (unit , file , status , operation , 
recno .buffer , size , ios ) 

character*(*)  f ile , status , operation 
character*20  tempop 
integer  unit , reel , ios , recno , size 
integer*2  buffer(siz2) 

tempop  -  operation 
call  casefold  (tempop)  1 

if  ( tempop, eq. 'OPEN' )  then 
if  (recl.gt.O)  then 
reel  -  size  *  2 

open  (unit-unit , file-file , status-status , access- ' direct ' , 
reel-reel , iostat-ios) 

endif 

else  if  ( tempop . eq . 'READ ' )  then 

read  (unit-unit , rec-recno , iostat-ios)  buffer 
else  if  ( tempop . eq . 'WRITE ' )  then 

write  (unit-unit , rec-recno, iostat-ios)  buffer 
endif 
return 
end 
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